\documentclass[]{report}

\usepackage{amsmath}
% uncomment if generating pdf documents
% comment if generating dvi documents
%\usepackage[pdftex]{color,graphicx}
%\usepackage[pdftex]{hyperref}
%\pdfpagewidth=\paperwidth
%\pdfpageheight=\paperheight

% comment if generating pdf documents
% uncomment if generating dvi documents
\usepackage[dvips]{graphicx}

% comment if using miktex(windows) or generating final release
% uncomment if using tetex(linux) and want to enable inverse search
%\usepackage[active]{srcltx}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\frontmatter
%\pagestyle{headings}
\title{Non-Topographic Photogrammetry}
\author{Richard G. Edwards}
\maketitle
\tableofcontents
\listoffigures
\listoftables


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\mainmatter
\chapter{The Photogrammetry Problem}
\section{Introduction}

The problem addressed in this paper is that of reconstructing a 3-d scene of object points from two or more photographs taken at different locations and orientations.  In order to solve such a problem a mathematical model that describes the physics will be derived.  As in any model there will be assumptions and simplifications made that will affect the overall accuracy of the solution.  The end goal is to take two or more photographs of an object and then find points on the surface of the object in the images, and then as output obtain where these points are located in 3-d space.  The camera parameters (location, orientation, etc.) will also be found in this process.
\section{Formulation of Local Mathematical Model}

Figure \ref{fig:collinear} shows an illustration of the photogrammetry model.  In the model, a physical object is represented as points which lie on the surface of the object, which will be known as the \emph{object points}.  The object is assumed to scatter all light from the surface and that the surface is not transparent, or in other words refraction in not considered in the model.  The model has one or more cameras that take two or more images of the object.  The camera is located at a displacement from the object with a certain orientation.  The camera is modeled as an optical system where all scattered light from the object passes through an image plane and converges to a single point, known as the \emph{perspective center}.  The \emph{image points} are the points on the image plane that correlate with the object points through the imaging process.

\begin{figure}
\begin{center}
\includegraphics[angle=0, width=1.0\textwidth]{collinear}
\end{center}
\caption{Mathematical model of imaging process with collinear rays.}\label{fig:collinear}
\end{figure}

A necessary step in solving such a problem is to describe the object points and camera parameters in a coordinate system.  The coordinate system is defined as a right handed Cartesian (rectangular) coordinate system and will be called the \emph{global} or \emph{object} coordinate system.  The coordinate system has three orthogonal basis vectors of unit length, which are represented as a symbol with a hat, such as $\hat{r}$. Orthogonality of basis vectors means that they have the property that the inner (dot) product of two orthogonal vectors is equal to zero.  Also the inner product of a basis vector with itself is equal to unity.  These properties are what defines the inner product operator.  The location of points in space are represented as an expansion of basis vectors with constant coefficients known as the the vector components.  The representation of a vector in this text will be a symbol with an arrow over the top of it, such as $\vec{r}$.  The arrow is a reminder that a vector has both a magnitude and a direction.  The camera interior parameters will be described in a local coordinate system called the \emph{image} coordinate system. In order to perform any mathematical operation on two or more vectors they must first be described in the same coordinate system.  This is accomplished by means of a transformation matrix.  

\subsection{Transformation Matrix}
A vector is used to represent a physical quantity that has both a magnitude and direction, such as displacement. It is a physical quantity that is the same regardless of the coordinate system used to describe it.  It can be expanded in any number of coordinate systems.  For example, the vector $\vec{v}$ is described in two different coordinate systems as follows.

\begin{align}
	\vec{v} = v_x\hat{x}+v_y\hat{y}+v_z\hat{z} = v_x'\hat{x}'+v_y'\hat{y}'+v_z'\hat{z}'\label{eqn:vectorexp}
\end{align}
It is worth reiterating that the vector represents the physics of the situation and does not change, but it can be described in any number of coordinate systems.  The inner product operator can be used to represent the primed vector components in terms of the unprimed basis vectors, thus expanding it in the other coordinate system.  The symbol used to represent the inner product will be a dot, for example $\hat{x}\cdot\hat{y}$.  This is accomplished as follows,
\begin{align}
v_x' &= \vec{v}\cdot\hat{x}'=(v_x\hat{x}+v_y\hat{y}+v_z\hat{z})\cdot\hat{x}'=(\hat{x}'\cdot\hat{x})v_x+(\hat{x}'\cdot\hat{y})v_y+(\hat{x}'\cdot\hat{z})v_z\\
v_y' &= \vec{v}\cdot\hat{y}'=(v_x\hat{x}+v_y\hat{y}+v_z\hat{z})\cdot\hat{y}'=(\hat{y}'\cdot\hat{x})v_x+(\hat{y}'\cdot\hat{y})v_y+(\hat{y}'\cdot\hat{z})v_z\\
v_z' &= \vec{v}\cdot\hat{z}'=(v_x\hat{x}+v_y\hat{y}+v_z\hat{z})\cdot\hat{z}'=(\hat{z}'\cdot\hat{x})v_x+(\hat{z}'\cdot\hat{y})v_y+(\hat{z}'\cdot\hat{z})v_z
\end{align}
This can be represented in matrix form as follows,
\begin{align}
\left[ \begin{array}{c} v_x' \\ v_y' \\ v_z' \end{array} \right] &=
\left[ \begin{array}{ccc}
\hat{x}'\cdot\hat{x} & \hat{x}'\cdot\hat{y} & \hat{x}'\cdot\hat{z} \\
\hat{y}'\cdot\hat{x} & \hat{y}'\cdot\hat{y} & \hat{y}'\cdot\hat{z} \\
\hat{z}'\cdot\hat{x} & \hat{z}'\cdot\hat{y} & \hat{z}'\cdot\hat{z}
\end{array} \right]
\left[ \begin{array}{c} v_x \\ v_y \\ v_z \end{array} \right]\label{eqn:matrixM}
\end{align}
or in more compact matrix notation as,
\begin{align}
\vec{v}'&= \mathbf{M}\vec{v}
\end{align}

The same procedure can be used to describe the unprimed vector components in terms of the primed basis vectors.  This leads to the following result,
\begin{align}
\left[ \begin{array}{c} v_x \\ v_y \\ v_z \end{array} \right] &=
\left[ \begin{array}{ccc}
\hat{x}\cdot\hat{x}' & \hat{x}\cdot\hat{y}' & \hat{x}\cdot\hat{z}' \\
\hat{y}\cdot\hat{x}' & \hat{y}\cdot\hat{y}' & \hat{y}\cdot\hat{z}' \\
\hat{z}\cdot\hat{x}' & \hat{z}\cdot\hat{y}' & \hat{z}\cdot\hat{z}'
\end{array} \right]
\left[ \begin{array}{c} v_x' \\ v_y' \\ v_z' \end{array} \right]\label{eqn:invM}
\end{align}
which when compared with (\ref{eqn:matrixM}) can be expressed as,
\begin{align}
\vec{v}&= \mathbf{M^T}\vec{v}'
\end{align}
where $^T$ is the transpose operator. Equations (\ref{eqn:matrixM}) and (\ref{eqn:invM}) are seen to be the inverse of each other which leads to the following result,
\begin{align}
\mathbf{M^T}=\mathbf{M^{-1}}
\end{align}
which shows that $\mathbf{M}$ is unitary.

The inner product of any two unit vectors is defined to be equal to the cosine of the angle between them.  
\begin{align}
	\hat{u}'\cdot\hat{u} = \cos(\theta_{u',u})
\end{align}
where $u=x,y,z$. This definition satisfies the conditions of orthognality. The matrix $\mathbf{M}$ can be cast into an alternate form more compatible with matrix notation as,
\begin{align}
\mathbf{M} = 
\left[ \begin{array}{ccc}
c_{1'1} & c_{1'2} & c_{1'3} \\
c_{2'1} & c_{2'2} & c_{2'3} \\
c_{3'1} & c_{3'2} & c_{3'3}
\end{array} \right]\label{eqn:M}
\end{align}
where $\{x,y,z\}\Leftrightarrow\{1,2,3\}$.  The way to interpret this notation is best described by example.  The expression of $c_{1'2}$ would be interpreted as the cosine of the angle between the unit vector $\hat{x}'$ and the unit vector $\hat{y}$.  Another example is that the expression of $c_{3'3}$ is interpreted as the cosine of the angle between the unit vector $\hat{z}'$ and the unit vector $\hat{z}$.

It is important to note that even though there are nine terms in the matrix there are only three independent variables that determine the matrix.  In general these variables will be $\alpha_x$, $\alpha_y$ and $\alpha_z$.  Depending on how the rotation matrix $\mathbf{M}$ is defined these independent variables will change.  There are many ways to define the rotation matrix, two of the more popular parameterizations are by means of \emph{Euler angles} and \emph{quaternion} rotations.

\subsubsection{Euler Angles}
Euler angles are parameterizations which include rotations around the three independent axes in a specified order.  Each of these rotations is illustrated in figure(\ref{fig:axis_rotations}).
\begin{figure}
\begin{center}
\includegraphics[angle=0, width=1.0\textwidth]{axis_rotations}
\end{center}
\caption{Euler rotations around (a) x-axis (b) y-axis and (c) z-axis.}\label{fig:axis_rotations}
\end{figure}
%\subsubsection{Rotation around x-axis}
Consider if the primed coordinate system is rotated around the x-axis by an angle of $\theta_x$.  By inspection (figure \ref{fig:axis_rotations} a), the matrix $\mathbf{M}$ from (\ref{eqn:M}) can be written as,
\begin{align}
\mathbf{M_x} = 
\left[ \begin{array}{ccc}
1 & 0 & 0 \\
0 & \cos(\theta_x) & \cos(\frac{\pi}{2}-\theta_x) \\
0 & \cos(\frac{\pi}{2}+\theta_x) & \cos(\theta_x)
\end{array} \right]
=\left[\begin{array}{ccc} 1&0&0\\0&\cos(\theta_x)&\sin(\theta_x)\\0&-\sin(\theta_x)&\cos(\theta_x) \end{array}\right]\label{eqn:Mx}
\end{align}
%\subsubsection{Rotation around y-axis}
By following a similar procedure a matrix describing a rotation around the y-axis is found to be,
\begin{align}
\mathbf{M_y} = 
\left[ \begin{array}{ccc}
\cos(\theta_y) & 0 & \cos(\frac{\pi}{2}+\theta_y) \\
0 & 1 & 0 \\
\cos(\frac{\pi}{2}-\theta_y) & 0 & \cos(\theta_y)
\end{array} \right]
=\left[\begin{array}{ccc} \cos(\theta_y)&0&-\sin(\theta_y)\\ 0&1&0\\ \sin(\theta_y)&0&\cos(\theta_y) \end{array}\right]\label{eqn:My}
\end{align}
%\subsubsection{Rotation around z-axis}
Finally, a matrix describing the rotation around the z-axis is,
\begin{align}
\mathbf{M_z} = 
\left[ \begin{array}{ccc}
\cos(\theta_z) & \cos(\frac{\pi}{2}-\theta_z) & 0 \\
\cos(\frac{\pi}{2}+\theta_z) & \cos(\theta_z) & 0\\
0 & 0 & 1
\end{array} \right]
=\left[\begin{array}{ccc} \cos(\theta_z)&\sin(\theta_z)&0\\-\sin(\theta_z)&\cos(\theta_z)&0 \\ 0&0&1 \end{array}\right]\label{eqn:Mz}
\end{align}
The overall orientation can be found by multiplying these three matrices in any arbitrary order.  The order that will be used is the following,
\begin{align}
\mathbf{M} = \mathbf{M_z}\mathbf{M_y}\mathbf{M_x}\label{eqn:MxMyMz}
\end{align}
Multiplying the matrices (\ref{eqn:Mx})--(\ref{eqn:Mz}) according to the order in (\ref{eqn:MxMyMz}) yields,
\begin{align}
\mathbf{M} = 
\left[ \begin{array}{ccc}
c_zc_y & s_xs_yc_z+c_xs_z & -c_xs_yc_z+s_xs_z \\
-c_ys_z & -s_xs_ys_z+c_xc_z & c_xs_ys_z+s_xc_z \\
s_y & -s_xc_y & c_xc_y
\end{array} \right]\label{eqn:orientationmatrix}
\end{align}
where $c_u=\cos(\theta_u)$, and  $s_u=\sin(\theta_u)$ and where $u = x,y,z$.

\subsubsection{Quaternion Rotations}
A \emph{quaternion} is defined by $q=q_0+q_1{i}+q_2{j}+q_3{k}$ where $q_0$, $q_1$, $q_2$ and $q_3$ are real numbers.  The \emph{addition} and \emph{subtraction} of two quaternions $a$ and $b$ is defined as,
\begin{align}
{a}\pm{b}&=(a_0+a_1i+a_2j+a_3k) \pm (b_0+b_1i+b_2j+b_3k) \\
&= (a_0\pm{b_0})+(a_1\pm{b_1})i+(a_2\pm{b_2})j+(a_3\pm{b_3})k\\
&=\sum_{n=0}^3{a_n\pm{b_n}}
\end{align}
Multiplication for primitive elements $i$, $j$ and $k$ is defined by $i^2=j^2=k^2=ijk=-1$, $ij=-ji=k$, $jk=-kj=i$, and $ki=-ik=j$. \emph{Multiplication} is defined as,
\begin{align}
ab &= (a_0+a_1i+a_2j+a_3k)(b_0+b_1i+b_2j+b_3k) \\
&= (a_0b_0-a_1b_1-a_2b_2-a_3b_3)\\
&= (a_0b_1+a_1b_0+a_2b_3-a_3b_2)i\\
&= (a_0b_2-a_1b_3+a_2b_0+a_3b_1)j\\
&= (a_0b_3+a_1b_2-a_2b_1+a_3b_0)k
\end{align}
Multiplication is not necessarily commutative $ab\ne{ba}$ but it is always associative $abc=(ab)c=a(bc)$.
The \emph{conjugate} of a quaternion is defined as,
\begin{align}
q^* &= (q_0+q_1i+q_2j+q_3k)^* = q_0-q_1i-q_2j-q_3k
\end{align}
which also satisfies the following properties $(q^*)^*=q$ and $(pq)^*=p^*q^*$. 
The \emph{norm} of a quaternion is defined as,
\begin{align}
N(q) = N(q_0+q_1i+q_2j+q_3k) = q_0^2+q_1^2+q_2^2+q_3^2
\end{align}
The norm is a real valued function which satisfies the following properties $N(q^*)=N(q)$ and $N(qp)=N(q)N(p)$. The \emph{inverse} of a quaternion is denoted $q^{-1}$ and has the property of $qq^{-1}=q^{-1}q=1$ and is defined as
\begin{align}
q^{-1} = \frac{q^*}{N(q)}
\end{align}
where division of a quaternion with a real number is defined as each component of the quaternion divided by the real number.  Another useful property is,
\begin{align}
q_0 = \frac{q+q^*}{2}
\end{align}
Quaternions can also be viewed as a scalar with a vector component $q=q_0+\vec{q}$ where $\vec{q} = q_1i+q_2j+q_3k$.  Quaternion multiplication can be defined in terms of vector multiplication as
\begin{align}
ab = (a_0b_0-\vec{a}\cdot\vec{b})+a_0\vec{b}+b_0\vec{a}+\vec{a}\times\vec{b}
\end{align}
If $\vec{a}\times\vec{b}=0$ (vectors are parallel) then the multiplication is commutative $ab=ba$.
A \emph{unit quaternion} is defined as a quaternion in which $N(q)=1$.  The inverse of unit quaternion and the multiplication of a unit quaternion is also a unit quaternion.  A unit quaternion can be defined as,
\begin{align}
q &= \cos(\psi)+\hat{u}\sin(\psi)
\end{align}
where $\hat{u}$ is a vector of unit length.  The unit vector can described in spherical coordinates  as $\hat{u} = \sin(\theta)\cos(\phi)i+\sin(\theta)\sin(\phi)j+cos(\theta)k$.  A unit quaternion can then be expressed in terms of the three indpendent variables $\psi$, $\theta$ and $\phi$ as,
\begin{align}
q &= c_{\psi} + s_{\psi}s_{\theta}c_{\phi}i + s_{\psi}s_{\theta}s_{\phi}j + s_{\psi}c_{\theta}k\\
[q_0, q_1, q_2, q_3] &= [c_{\psi}, s_{\psi}s_{\theta}c_{\phi}, s_{\psi}s_{\theta}s_{\phi}, s_{\psi}c_{\theta}]\label{eqn:unitquaternion}
\end{align}
This quaternion is interpreted as a rotation about the unit vector $\hat{u}$ by an angle $\omega=2\psi$.
A vector $\vec{v}$ that has been rotated by a unit quaternion is given by,
\begin{align}
v' = qvq^*
\end{align}
where $v=v_0+\vec{v}$.
Conversion from a quaternion to a rotation matrix is given as,
\begin{align}
\mathbf{M} = 
\left
[ \begin{array}{ccc}
	q_0^2+q_1^2-q_2^2-q_3^2 & 2q_1q_2-2q_0q_3 & 2q_1q_3+2q_0q_2\\
	2q_1q_2+2q_0q_3 & q_0^2-q_1^2+q_2^2-q_3^2 & 2q_2q_3-2q_0q_1\\
	2q_1q_3-2q_0q_2 & 2q_2q_3+2q_0q_1 & q_0^2-q_1^2-q_2^2+q_3^2
\end{array}
\right]\label{eqn:Mquaternion}
\end{align}

\subsection{Collinearity Equation}

The mathematical expression that the object point, the image point and the perspective center of the camera ideally lie in a straight line is known as the \emph{collinearity equation}. Referring to figure \ref{fig:collinear}, the vector $\vec{u}-\vec{v}$ is collinear with the vector $\vec{r}-\vec{t}$ if they have the same direction.  In order for this to be true, only the direction of the vectors need to be equal, where the  magnitudes are proportional by an arbitrary scaling constant $k_s$.  Since each vector is described in a different coordinate system, one needs to be rotated into the coordinate system of the other by means of a rotation matrix.  This concept expressed mathematically in vector form is,
\begin{align}
\vec{u}-\vec{v} = k_s\mathbf{M}(\vec{r}-\vec{t})\label{eqn:collinear_vectorform}
\end{align}
Expanding (\ref{eqn:collinear_vectorform}) yields three equations,
\begin{align}
u_x-v_x &= k_s[m_{11}(r_x-t_x)+m_{12}(r_y-t_y)+m_{13}(r_z-t_z)]\label{eqn:collinear_expanded_x}\\
u_y-v_y &= k_s[m_{21}(r_x-t_x)+m_{22}(r_y-t_y)+m_{23}(r_z-t_z)]\label{eqn:collinear_expanded_y}\\
   -v_z &= k_s[m_{31}(r_x-t_x)+m_{32}(r_y-t_y)+m_{33}(r_z-t_z)]\label{eqn:collinear_expanded_z}
\end{align}
The scaling constant $k_s$ can be eliminated by solving for it in (\ref{eqn:collinear_expanded_z}) to get,
\begin{align}
k_s = -v_z[m_{31}(r_x-t_x)+m_{32}(r_y-t_y)+m_{33}(r_z-t_z)]^{-1}
\end{align}
and then substituting into (\ref{eqn:collinear_expanded_x}) and (\ref{eqn:collinear_expanded_y}) yields the collinearity equations in functional form,
\begin{align}
f_a(p) &= u_x-v_x+v_z\frac{m_{11}(r_x-t_x)+m_{12}(r_y-t_y)+m_{13}(r_z-t_z)}{m_{31}(r_x-t_x)+m_{32}(r_y-t_y)+m_{33}(r_z-t_z)}\label{eqn:collinear1}\\
f_b(p) &= u_y-v_y+v_z\frac{m_{21}(r_x-t_x)+m_{22}(r_y-t_y)+m_{23}(r_z-t_z)}{m_{31}(r_x-t_x)+m_{32}(r_y-t_y)+m_{33}(r_z-t_z)}\label{eqn:collinear2} 
\end{align}
where,
\begin{align}
p=\{r_x,r_y,r_z,t_x,t_y,t_z,v_x,v_y,v_z,\alpha_x,\alpha_y,\alpha_z,u_x,u_y\}\label{eqn:photogramset}
\end{align}
These fourteen independent variables are separated into subsets of,
\begin{itemize}
\item object parameters ${r_x, r_y, r_z}$
\item camera parameters ${v_x, v_y, v_z}$
\item image parameters ${t_x, t_y, t_z, \alpha_x, \alpha_y, \alpha_z, u_x, u_y}$
\end{itemize} 
The reason for separating the camera parameters from the image parameters is that different images can be taken with the same camera. The descriptions of each variable along with the coordinate system in which it is expanded is found in table \ref{tab:DescriptionOfIndependentVariables}.
\begin{table}
	\caption{Description of Independent Variables}
	\centering
	\begin{tabular}{|l|l|l|}
	\hline
	Variables  & Description & Coordinate system\\
	\hline
	\hline
	$r_x$, $r_y$, $r_z$ & location of object points & object\\
	\hline
	$t_x$, $t_y$, $t_z$ & location of image & object\\
	\hline
	$u_x$, $u_y$ & location of image points & image\\
	\hline
	$v_x$, $v_y$, $v_z$   & principal center & image\\
	\hline
	$\alpha_x$, $\alpha_y$, $\alpha_z$ & rotation matrix parameters & image\\
	\hline
	\end{tabular}
	\label{tab:DescriptionOfIndependentVariables}
\end{table}

There are many ways of combining (\ref{eqn:collinear1}) and (\ref{eqn:collinear2}) into a single function.  One common way to combine them is to take the root of the mean of the square (RMS) of each function,
\begin{align}
f(p) = \sqrt{\frac{{f_a(p)}^2+{f_b(p)}^2}{2}}\label{eqn:collinearcombine}
\end{align}
The reason for combining the two equations into one single equation is that both of them need to be satisfied, in which they should ideally be equal to zero.  Because of approximations in the model they will never be equal to zero but equal to a small error.  The global error from both functions will be minimized to obtain the best overall solution. Therefore, the error of the individual functions need to be combined in some way to show how they relate to each other.

\subsection{Gradient of Collinearity Equations}
The gradient of a scalar function is defined as,
\begin{align}
\nabla{f(x)}=\sum_n\frac{\partial{f}}{\partial{x_i}}\hat{x}_i
\end{align}
where $p$ is defined in (\ref{eqn:photogramset}).  The gradient of the scalar functions (\ref{eqn:collinear1}) and (\ref{eqn:collinear2}) can be found analytically. In order to accomplish this the following auxiliary equation is defined,
\begin{align}
\vec{w} = \mathbf{M}(\vec{r}-\vec{t})\label{eqn:vecw}
\end{align}
which when expanded yields,
\begin{align}
\left[ \begin{array}{c}
w_x\\w_y\\w_z
\end{array} \right]
=\left[ \begin{array}{ccc}
m_{11}(r_x-t_x)+m_{12}(r_y-t_y)+m_{13}(r_z-t_z)\\
m_{21}(r_x-t_x)+m_{22}(r_y-t_y)+m_{23}(r_z-t_z)\\
m_{31}(r_x-t_x)+m_{32}(r_y-t_y)+m_{33}(r_z-t_z)
\end{array} \right]
\end{align}
Substituting into (\ref{eqn:collinear1}) and (\ref{eqn:collinear2}) gives,
\begin{align}
f_a(p) &= u_x-v_x+v_z\frac{w_x}{w_z}\label{eqn:collineara_1}\\
f_b(p) &= u_y-v_y+v_z\frac{w_y}{w_z}\label{eqn:collinearb_1}
\end{align}
By making use of the chain rule, (\ref{eqn:collineara_1}) and ({\ref{eqn:collinearb_1}) are differentiated with respect to an arbitrary variable $p_n$,
\begin{align}
\frac{\partial{f_a}}{\partial{p_n}}=\frac{\partial{u_x}}{\partial{p_n}}-\frac{\partial{v_x}}{\partial{p_n}}+\frac{w_x}{w_z}\frac{\partial{v_z}}{\partial{p_n}}+\frac{v_z}{w_z}\left(\frac{\partial{w_x}}{\partial{p_n}}-\frac{w_x}{w_z}\frac{\partial{w_z}}{\partial{p_n}}\right)\label{eqn:dcollineara}\\
\frac{\partial{f_b}}{\partial{p_n}}=\frac{\partial{u_y}}{\partial{p_n}}-\frac{\partial{v_y}}{\partial{p_n}}+\frac{w_y}{w_z}\frac{\partial{v_z}}{\partial{p_n}}+\frac{v_z}{w_z}\left(\frac{\partial{w_y}}{\partial{p_n}}-\frac{w_y}{w_z}\frac{\partial{w_z}}{\partial{p_n}}\right)\label{eqn:dcollinearb}
\end{align}
where $p_n$ is equal to any of the variables in (\ref{eqn:photogramset}).
Using (\ref{eqn:dcollineara}) and (\ref{eqn:dcollinearb}) and letting $p_n=\{u_x, u_y\}$ the partial derivatives of the collinearity equations with respect to image locations $\vec{u}$ are found to be,
\begin{align}
\left[ \begin{array}{cc}
\frac{\partial{f_a}}{\partial{u_x}} & \frac{\partial{f_a}}{\partial{u_y}} \\
\frac{\partial{f_b}}{\partial{u_x}} & \frac{\partial{f_b}}{\partial{u_y}}
\end{array} \right]
&=\left[ \begin{array}{cc}
1 & 0 \\
0 & 1
\end{array} \right]
\end{align}
Next, with respect to the principal point location $\vec{v}$,
\begin{align}
\left[ \begin{array}{ccc}
\frac{\partial{f_a}}{\partial{v_x}} & \frac{\partial{f_a}}{\partial{v_y}} & \frac{\partial{f_a}}{\partial{v_z}}\\
\frac{\partial{f_b}}{\partial{v_x}} & \frac{\partial{f_b}}{\partial{v_y}} & \frac{\partial{f_a}}{\partial{v_z}}
\end{array} \right]
=\left[ \begin{array}{ccc}
-1 & 0 & {w_x}/{w_z} \\
0 & -1 & {w_y}/{w_z}
\end{array} \right]
\end{align}
Then with respect to the object point $\vec{r}$ and camera position $\vec{t}$,
\begin{align}
\left[\begin{array}{ccc}
\frac{\partial{f_a}}{\partial{r_x}}\\
\frac{\partial{f_a}}{\partial{r_y}}\\
\frac{\partial{f_a}}{\partial{r_z}}
\end{array}\right]&=-
\left[\begin{array}{ccc}
\frac{\partial{f_a}}{\partial{t_x}}\\
\frac{\partial{f_a}}{\partial{t_y}}\\
\frac{\partial{f_a}}{\partial{t_z}}
\end{array}\right]
=\frac{v_z}{w_z}
\left[\begin{array}{ccc}
m_{11}-m_{31}({w_x}/{w_z})\\
m_{12}-m_{32}({w_x}/{w_z})\\
m_{13}-m_{33}({w_x}/{w_z})
\end{array}\right]\\
\left[\begin{array}{ccc}
\frac{\partial{f_b}}{\partial{r_x}}\\
\frac{\partial{f_b}}{\partial{r_y}}\\
\frac{\partial{f_b}}{\partial{r_z}}
\end{array}\right]&=-
\left[\begin{array}{ccc}
\frac{\partial{f_b}}{\partial{t_x}}\\
\frac{\partial{f_b}}{\partial{t_y}}\\
\frac{\partial{f_b}}{\partial{t_z}}
\end{array}\right]
=\frac{v_z}{w_z}
\left[\begin{array}{ccc}
m_{21}-m_{31}({w_y}/{w_z})\\
m_{22}-m_{32}({w_y}/{w_z})\\
m_{23}-m_{33}({w_y}/{w_z})
\end{array}\right]
\end{align}
The partial derivatives with respect to the image orientation {$\alpha_x$, $\alpha_y$, $\alpha_z$} are calculated by using (\ref{eqn:dcollineara}) and (\ref{eqn:dcollinearb}) with $p_n=\{\alpha_x,\alpha_y,\alpha_z\}$ to get,
\begin{align}
\frac{\partial{f_a}}{\partial{\alpha_n}}=\frac{v_z}{w_z}\left(\frac{\partial{w_x}}{\partial{\alpha_n}}-\frac{w_x}{w_z}\frac{\partial{w_z}}{\partial{\alpha_n}}\right)\label{eqn:dcollineara_alpha}\\
\frac{\partial{f_b}}{\partial{\alpha_n}}=\frac{v_z}{w_z}\left(\frac{\partial{w_y}}{\partial{\alpha_n}}-\frac{w_y}{w_z}\frac{\partial{w_z}}{\partial{\alpha_n}}\right)\label{eqn:dcollinearb_alpha}
\end{align}
with $n=x,y,z$. The partial derivatives of $\vec{w}$ with respect to the image orientation is found by differentiating (\ref{eqn:vecw}) to get, 
\begin{align}
\frac{\partial\vec{w}}{\partial\alpha_n} = \frac{\partial\mathbf{M}}{\partial\alpha_n}(\vec{r}-\vec{t})\label{eqn:dvecw}
\end{align}
with $n=x,y,z$ and where $\vec{r}-\vec{t}$ is a constant. It should be noted that the above results are general with respect to the orientation matrix.  As stated earlier there are numerous ways of specifying the orientation matrix.  For example the order in which they are multiplied could be changed or quaternions could be used. From (\ref{eqn:dvecw}) it can be seen that the only remaining thing to find is the partial derivatives of the orientation matrix $\mathbf{M}$ with respect to the orientation variables $\alpha_x$, $\alpha_y$ and $\alpha_z$. 

\subsubsection{Differentiation of rotation matrix with respect to Euler angles}
When using Euler angles the independent variables are angular rotations about the x-axis ($\theta_x$), the y-axis ($\theta_y$) and the z-axis ($\theta_z$) according to the order specified in (\ref{eqn:MxMyMz}).
Differentiating (\ref{eqn:orientationmatrix}) yields,
\begin{align}
\frac{\partial\mathbf{M}}{\partial{\theta_x}}&=
\left[ \begin{array}{ccc}
0 & c_xs_yc_z-s_xs_z & s_xs_yc_z+c_xs_z \\
0 & -c_xs_ys_z-s_xc_z & -s_xs_ys_z+c_xc_z \\
0 & -c_xc_y & -s_xc_y
\end{array} \right]\label{eqn:dM_dthetax}\\
&=\left[ \begin{array}{ccc}
0 & -m_{13} & m_{12}\\
0 & -m_{23} & m_{22}\\
0 & -m_{33} & m_{32}
\end{array} \right]\label{eqn:dM_dthetax2}
\end{align}
\begin{align}
\frac{\partial\mathbf{M}}{\partial{\theta_y}} &= 
\left[ \begin{array}{ccc}
-c_zs_y & s_xc_yc_z & -c_xc_yc_z \\
s_ys_z & -s_xc_ys_z & c_xc_ys_z \\
c_y & s_xs_y & -c_xs_y
\end{array} \right]\label{eqn:dM_dthetay}\\
&=\left[ \begin{array}{ccc}
-c_zm_{31} & -c_zm_{32} & -c_zm_{33}\\
s_zm_{31} & s_zm_{32} & s_zm_{33}\\
c_zm_{11}-s_zm_{21} & c_zm_{12}-s_zm_{22} &c_zm_{13}-s_zm_{23}
\end{array} \right]\label{eqn:dM_dthetay2}
\end{align}
\begin{align}
\frac{\partial\mathbf{M}}{\partial{\theta_z}} &= 
\left[ \begin{array}{ccc}
-s_zc_y & -s_xs_ys_z+c_xc_z & c_xs_ys_z+s_xc_z \\
-c_yc_z & -s_xs_yc_z-c_xs_z & c_xs_yc_z-s_xs_z \\
0 & 0 & 0
\end{array} \right]\label{eqn:dM_dthetaz}\\
&=\left[ \begin{array}{ccc}
m_{21} & m_{22} & m_{23}\\
-m_{11} & -m_{12} & -m_{13}\\
0 & 0 & 0
\end{array} \right]\label{eqn:dM_dthetaz2}
\end{align}
It should also be noted that since there are zeros in the partial derivatives of the orientation matrix, it is more computationally efficient to carry out the multiplications.  By doing this, the partial derivatives of $f_a$ and $f_b$ with respect to $\theta_x$, $\theta_y$ and $\theta_z$ are found to be,
\begin{align}
\begin{split}
\frac{\partial{f_a}}{\partial{\theta_x}}=&\frac{v_z}{w_z}\biggl[m_{12}(r_z-t_z)-m_{13}(r_y-t_y)\\
&-\frac{w_x}{w_z}[m_{32}(r_z-t_z)-m_{33}(r_y-t_y)]\biggr]
\end{split}\\
\frac{\partial{f_a}}{\partial{\theta_y}}=&\frac{v_z}{w_z}\left[-c_zw_z-\frac{w_x}{w_z}(c_zw_x-s_zw_y)\right]\\
\frac{\partial{f_a}}{\partial{\theta_z}}=&\frac{v_zw_y}{w_z}
\end{align}
\begin{align}
\begin{split}
\frac{\partial{f_b}}{\partial{\theta_x}}=&\frac{v_z}{w_z}\biggl[m_{22}(r_z-t_z)-m_{23}(r_y-t_y)\\
&-\frac{w_y}{w_z}[m_{32}(r_z-t_z)-m_{33}(r_y-t_y)]\biggr]
\end{split}\\
\frac{\partial{f_b}}{\partial{\theta_y}}=&\frac{v_z}{w_z}\left[s_zw_z-\frac{w_y}{w_z}(c_zw_x-s_zw_y)\right]\\
\frac{\partial{f_b}}{\partial{\theta_z}}=&\frac{-v_zw_x}{w_z}
\end{align}

\subsubsection{Differentiation of rotation matrix from quaternions}
The partial derivatives of the rotation matrix $\mathbf{M}$ derived from quaternions is calculated by taking the partial derivatives with respect the three independent variables $\psi$, $\theta$ and $\phi$ of (\ref{eqn:unitquaternion}).  First, with respect to $\psi$,
\begin{align}
\frac{\partial{q}}{\partial{\psi}} &= 
-s_{\psi} + c_{\psi}s_{\theta}c_{\phi}i + c_{\psi}s_{\theta}s_{\phi}j + c_{\psi}c_{\theta}k\\
[q_0',q_1',q_2',q_3']&=[-s_{\psi}, c_{\psi}s_{\theta}c_{\phi}, c_{\psi}s_{\theta}s_{\phi}, c_{\psi}c_{\theta}]
\end{align}
where $'=\frac{\partial}{\partial{\psi}}$.  Then with respect to $\theta$,
\begin{align}
\frac{\partial{q}}{\partial{\theta}} &= 
s_{\psi}c_{\theta}c_{\phi}i + s_{\psi}c_{\theta}s_{\phi}j - s_{\psi}s_{\theta}k\\
[q_0',q_1',q_2',q_3']&=[0,s_{\psi}c_{\theta}c_{\phi}, s_{\psi}c_{\theta}s_{\phi}, - s_{\psi}s_{\theta}]
\end{align}
where $'=\frac{\partial}{\partial{\theta}}$.  And finally with respect to $\phi$,
\begin{align}
\frac{\partial{q}}{\partial{\phi}} &= 
-s_{\psi}s_{\theta}s_{\phi}i + s_{\psi}s_{\theta}c_{\phi}j\\
[q_0',q_1',q_2',q_3']&=[0,-s_{\psi}s_{\theta}s_{\phi}, s_{\psi}s_{\theta}c_{\phi},0]
\end{align}
where $'=\frac{\partial}{\partial{\phi}}$.  The partial derivatives of the rotation matrix are then found by differentiating (\ref{eqn:Mquaternion}) to get,
\begin{align}
m_{11}' &= 2(q_0q_0'+q_1'q_1-q_2q_2'-q_3'q_3)\\ 
m_{12}' &= 2(q_1q_2'+q_1'q_2-q_0q_3'-q_0'q_3)\\
m_{13}' &= 2(q_3q_1'+q_3'q_1+q_0q_2'+q_0'q_2)\\
m_{21}' &= 2(q_1q_2'+q_1'q_2+q_0q_3'+q_0'q_3)\\
m_{22}' &= 2(q_0q_0'-q_1'q_1+q_2q_2'-q_3'q_3)\\
m_{23}' &= 2(q_2q_3'+q_2'q_3-q_0q_1'-q_0'q_1)\\
m_{31}' &= 2(q_3q_1'+q_3'q_1-q_0q_2'-q_0'q_2)\\
m_{32}' &= 2(q_2q_3'+q_2'q_3+q_0q_1'+q_0'q_1)\\
m_{33}' &= 2(q_0q_0'-q_1'q_1-q_2q_2'+q_3'q_3)
\end{align}
where $'=\frac{\partial}{\partial{\psi}},\frac{\partial}{\partial{\theta}},\frac{\partial}{\partial{\phi}}$.

\subsubsection{Differentiation of combined collinearity equations}

The partial derivatives of the combined collinearity equation of (\ref{eqn:collinearcombine}) is found to be,
\begin{align}
\frac{\partial{f}}{\partial{p_n}}=\frac{1}{2f}\left[f_a\frac{\partial{f_a}}{\partial{p_n}}+f_b\frac{\partial{f_b}}{\partial{p_n}}\right]\label{eqn:dcollinearcombine}
\end{align}

\subsection{Hessian matrix}
The \emph{Hessian matrix} of a scalar function is defined as,
\begin{align}
f''(p)&\equiv\sum_i\sum_j\frac{\partial^2f(p)}{\partial{x_i\partial{x_j}}}\hat{x_{ij}}=
\left[ \begin{array}{cccc}
\frac{\partial^2{f}}{\partial{p_1}\partial{p_1}} & \frac{\partial^2{f}}{\partial{p_1}\partial{p_2}} & \ldots & \frac{\partial^2{f}}{\partial{p_1}\partial{p_N}} \\
\frac{\partial^2{f}}{\partial{p_2}\partial{p_1}} & \frac{\partial^2{f}}{\partial{p_2}\partial{p_2}} & & \frac{\partial^2{f}}{\partial{p_2}\partial{p_N}} \\
\vdots &  & \ddots & \vdots \\
\frac{\partial^2{f}}{\partial{p_M}\partial{p_1}} & \frac{\partial^2{f}}{\partial{p_M}\partial{p_2}} & \ldots & \frac{\partial^2{f}}{\partial{p_M}\partial{p_N}} \\
\end{array} \right]\label{eqn:Hessian}
\end{align}
This Hessian matrix can be seen to be the Jacobian of the gradient.  Differentiating (\ref{eqn:dcollineara}) and (\ref{eqn:dcollinearb}) with respect to an arbitrary variable $q$ and using the chain rule yields,
\begin{align}
\frac{\partial^2{f_a}}{\partial{q}\partial{p}}&=\frac{w_z}{v_z}\left[\left(\frac{\partial}{\partial{q}}\frac{v_z}{w_z}\right)\frac{\partial{f_a}}{\partial{p}}+\left(\frac{\partial}{\partial{p}}\frac{v_z}{w_z}\right)\frac{\partial{f_a}}{\partial{q}}\right]+\frac{v_z}{w_z}\left[\frac{\partial^2{w_x}}{\partial{q}\partial{p}}-\frac{w_x}{w_z}\frac{\partial^2{w_z}}{\partial{q}\partial{p}}\right]\label{eqn:ddcollineara}\\
\frac{\partial^2{f_b}}{\partial{q}\partial{p}}&=\frac{w_z}{v_z}\left[\left(\frac{\partial}{\partial{q}}\frac{v_z}{w_z}\right)\frac{\partial{f_b}}{\partial{p}}+\left(\frac{\partial}{\partial{p}}\frac{v_z}{w_z}\right)\frac{\partial{f_b}}{\partial{q}}\right]+\frac{v_z}{w_z}\left[\frac{\partial^2{w_y}}{\partial{q}\partial{p}}-\frac{w_y}{w_z}\frac{\partial^2{w_z}}{\partial{q}\partial{p}}\right]\label{eqn:ddcollinearb}
\end{align}
where, 
\begin{align}
\frac{\partial}{\partial{p}}\frac{v_z}{w_z}&=\frac{1}{w_z}\left(\frac{\partial{v_z}}{\partial{p}}-\frac{v_z}{w_z}\frac{\partial{w_z}}{\partial{p}}\right)\label{eqn:dvz_wz_p}\\
\frac{\partial}{\partial{q}}\frac{v_z}{w_z}&=\frac{1}{w_z}\left(\frac{\partial{v_z}}{\partial{q}}-\frac{v_z}{w_z}\frac{\partial{w_z}}{\partial{q}}\right)\label{eqn:dvz_wz_q}\\
\end{align}
When $p=q$ then (\ref{eqn:ddcollineara}) and (\ref{eqn:ddcollinearb}) reduce to,
\begin{align}
\frac{\partial^2{f_a}}{\partial{q}\partial{p}}&=\frac{2w_z}{v_z}\left[\left(\frac{\partial}{\partial{p}}\frac{v_z}{w_z}\right)\frac{\partial{f_a}}{\partial{p}}\right]+\frac{v_z}{w_z}\left[\frac{\partial^2{w_x}}{\partial{p}^2}-\frac{w_x}{w_z}\frac{\partial^2{w_z}}{\partial{p}^2}\right]\label{eqn:ddcollinearapq}\\
\frac{\partial^2{f_b}}{\partial{q}\partial{p}}&=\frac{2w_z}{v_z}\left[\left(\frac{\partial}{\partial{p}}\frac{v_z}{w_z}\right)\frac{\partial{f_b}}{\partial{p}}\right]+\frac{v_z}{w_z}\left[\frac{\partial^2{w_y}}{\partial{p}^2}-\frac{w_y}{w_z}\frac{\partial^2{w_z}}{\partial{p}^2}\right]\label{eqn:ddcollinearbpq}
\end{align}
Differentiating $\vec{w}$ yields,
\begin{align}
\frac{\partial\vec{w}}{\partial{u_x}}=\frac{\partial\vec{w}}{\partial{u_y}}=\frac{\partial\vec{w}}{\partial{v_x}}=\frac{\partial\vec{w}}{\partial{v_y}}=\frac{\partial\vec{w}}{\partial{v_z}}=0
\end{align}
\begin{align}
	\frac{\partial\vec{w}}{\partial{r_x}}&=-\frac{\partial\vec{w}}{\partial{t_x}}=\textbf{M}\hat{x}\\
	\frac{\partial\vec{w}}{\partial{r_y}}&=-\frac{\partial\vec{w}}{\partial{t_y}}=\textbf{M}\hat{y}\\
	\frac{\partial\vec{w}}{\partial{r_z}}&=-\frac{\partial\vec{w}}{\partial{t_z}}=\textbf{M}\hat{z}
\end{align}
where $\hat{x}^T=[1,0,0]$,$\hat{y}^T=[0,1,0]$ and $\hat{z}^T=[0,0,1]$.
\begin{align}
	\frac{\partial\vec{w}}{\partial{q_w}}&=\frac{\partial\textbf{M}}{\partial{q_w}}(\vec{r}-\vec{t})\\
	\frac{\partial\vec{w}}{\partial{q_x}}&=\frac{\partial\textbf{M}}{\partial{q_x}}(\vec{r}-\vec{t})\\
	\frac{\partial\vec{w}}{\partial{q_y}}&=\frac{\partial\textbf{M}}{\partial{q_y}}(\vec{r}-\vec{t})\\
	\frac{\partial\vec{w}}{\partial{q_z}}&=\frac{\partial\textbf{M}}{\partial{q_z}}(\vec{r}-\vec{t})
\end{align}
where,
\begin{align}
\frac{\partial\mathbf{M}}{\partial{q_w}} = 2
\left
[ \begin{array}{ccc}
	 q_0 & -q_3 &  q_2\\
	 q_3 &  q_0 & -q_1\\
	-q_2 &  q_1 &  q_0
\end{array}
\right]\label{eqn:dMqw}
\end{align}
\begin{align}
\frac{\partial\mathbf{M}}{\partial{q_x}} = 2
\left
[ \begin{array}{ccc}
	q_1 &  q_2 &  q_3\\
	q_2 & -q_1 & -q_0\\
	q_3 &  q_0 & -q_1
\end{array}
\right]\label{eqn:dMqx}
\end{align}
\begin{align}
\frac{\partial\mathbf{M}}{\partial{q_y}} = 2
\left
[ \begin{array}{ccc}
	-q_2 & q_1 &  q_0\\
	 q_1 & q_2 &  q_3\\
	-q_0 & q_3 & -q_2
\end{array}
\right]\label{eqn:dMqy}
\end{align}
\begin{align}
\frac{\partial\mathbf{M}}{\partial{q_z}} = 2
\left
[ \begin{array}{ccc}
	-q_3 & -q_0 & q_1\\
	 q_0 & -q_3 & q_2\\
	 q_1 &  q_2 & q_3
\end{array}
\right]\label{eqn:dMqz}
\end{align}
Differentiating $\vec{w}$ again gives the following partial differential equations,
\begin{align}
\frac{\partial\vec{w}}{\partial{q_m}\partial{r_x}}&=-\frac{\partial\vec{w}}{\partial{q_m}\partial{t_x}}=\frac{\partial\mathbf{M}}{\partial{q_m}}\hat{x}\\
\frac{\partial\vec{w}}{\partial{q_m}\partial{r_y}}&=-\frac{\partial\vec{w}}{\partial{q_m}\partial{t_y}}=\frac{\partial\mathbf{M}}{\partial{q_m}}\hat{y}\\
\frac{\partial\vec{w}}{\partial{q_m}\partial{r_z}}&=-\frac{\partial\vec{w}}{\partial{q_m}\partial{t_z}}=\frac{\partial\mathbf{M}}{\partial{q_m}}\hat{z}
\end{align}
where $m=w,x,y,z$.
\begin{align}
\frac{\partial^2\mathbf{M}}{\partial{q_m}\partial{q_n}} &= 2
\left
[ \begin{array}{cccc}
	\frac{\partial^2\mathbf{M}}{\partial{q_w}\partial{q_w}} & \frac{\partial^2\mathbf{M}}{\partial{q_w}\partial{q_x}} & \frac{\partial^2\mathbf{M}}{\partial{q_w}\partial{q_y}} & \frac{\partial^2\mathbf{M}}{\partial{q_w}\partial{q_z}}\\
	\frac{\partial^2\mathbf{M}}{\partial{q_x}\partial{q_w}} & \frac{\partial^2\mathbf{M}}{\partial{q_x}\partial{q_x}} & \frac{\partial^2\mathbf{M}}{\partial{q_x}\partial{q_y}} & \frac{\partial^2\mathbf{M}}{\partial{q_x}\partial{q_z}}\\
	\frac{\partial^2\mathbf{M}}{\partial{q_y}\partial{q_w}} & \frac{\partial^2\mathbf{M}}{\partial{q_y}\partial{q_x}} & \frac{\partial^2\mathbf{M}}{\partial{q_y}\partial{q_y}} & \frac{\partial^2\mathbf{M}}{\partial{q_y}\partial{q_z}}\\
	\frac{\partial^2\mathbf{M}}{\partial{q_z}\partial{q_w}} & \frac{\partial^2\mathbf{M}}{\partial{q_z}\partial{q_x}} & \frac{\partial^2\mathbf{M}}{\partial{q_z}\partial{q_y}} & \frac{\partial^2\mathbf{M}}{\partial{q_z}\partial{q_z}}
\end{array}
\right]\\
&=2\left[ \begin{array}{ccccccccccccccc}
	  1 &  0 &  0 & & 0 &  0 &  0 & &  0 & 0 &  1 & &  0 & -1 & 0\\
	  0 &  1 &  0 & & 0 &  0 & -1 & &  0 & 0 &  0 & &  1 &  0 & 0\\
	  0 &  0 &  1 & & 0 &  1 &  0 & & -1 & 0 &  0 & &  0 &  0 & 0\\
	  \\
	  0 &  0 &  0 & & 1 &  0 &  0 & &  0 & 1 &  0 & &  0 &  0 & 1\\
	  0 &  0 & -1 & & 0 & -1 &  0 & &  1 & 0 &  0 & &  0 &  0 & 0\\
	  0 &  1 &  0 & & 0 &  0 & -1 & &  0 & 0 &  0 & &  1 &  0 & 0\\
	  \\
	  0 &  0 &  1 & & 0 &  1 &  0 & & -1 & 0 &  0 & &  0 &  0 & 0\\
	  0 &  0 &  0 & & 1 &  0 &  0 & &  0 & 1 &  0 & &  0 &  0 & 1\\
	 -1 &  0 &  0 & & 0 &  0 &  0 & &  0 & 0 & -1 & &  0 &  1 & 0\\
	 \\
	  0 & -1 &  0 & & 0 &  0 &  1 & &  0 & 0 &  0 & & -1 &  0 & 0\\
	  1 &  0 &  0 & & 0 &  0 &  0 & &  0 & 0 &  1 & &  0 & -1 & 0\\
	  0 &  0 &  0 & & 1 &  0 &  0 & &  0 & 1 &  0 & &  0 &  0 & 1\\
\end{array}\right]
\end{align}
When $q$ or $p=u_x, u_y, v_x, v_y$ then (\ref{eqn:ddcollineara}) and (\ref{eqn:ddcollinearb}) equate to,
\begin{align}
\frac{\partial^2{f_a}}{\partial{q}\partial{p}}=\frac{\partial^2{f_b}}{\partial{q}\partial{p}}=0
\end{align}
Differentiating with respect to the variable $v_z$, equations (\ref{eqn:ddcollineara}) and (\ref{eqn:ddcollinearb}) reduce to,
\begin{align}
\frac{\partial^2{f_a}}{\partial{q}\partial{v_z}}&=\frac{1}{v_z}\frac{\partial{f_a}}{\partial{q}}\\
\frac{\partial^2{f_b}}{\partial{q}\partial{v_z}}&=\frac{1}{v_z}\frac{\partial{f_b}}{\partial{q}}\\
\frac{\partial^2{f_a}}{\partial{v_z}^2}&=\frac{\partial^2{f_b}}{\partial{v_z}^2}=0
\end{align}
where $q={r_x, r_y, r_z, t_x, t_y, t_z, q_w, q_x, q_y, q_z}$.  Letting $q$ or $p=q_w, q_x, q_y, q_z$,

\begin{align}
\frac{\partial^2{f_a}}{\partial{q_w^2}}&=\frac{-2}{w_z}\frac{\partial{f_a}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_w}}+\frac{2v_z}{w_z}\left[(r_x-t_x)-\frac{w_x}{w_z}(r_z-t_z)\right]\label{eqn:df_aq_wq_w}\\
\frac{\partial^2{f_b}}{\partial{q_w^2}}&=\frac{-2}{w_z}\frac{\partial{f_b}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_w}}+\frac{2v_z}{w_z}\left[(r_y-t_y)-\frac{w_y}{w_z}(r_z-t_z)\right]\label{eqn:df_bq_wq_w}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_x^2}}&=\frac{-2}{w_z}\frac{\partial{f_a}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_x}}+\frac{2v_z}{w_z}\left[(r_x+t_x)+\frac{w_x}{w_z}(r_z-t_z)\right]\label{eqn:df_aq_xq_x}\\
\frac{\partial^2{f_b}}{\partial{q_x^2}}&=\frac{-2}{w_z}\frac{\partial{f_b}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_x}}+\frac{2v_z}{w_z}\left[-(r_y+t_y)+\frac{w_y}{w_z}(r_z-t_z)\right]\label{eqn:df_bq_xq_x}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_y^2}}&=\frac{-2}{w_z}\frac{\partial{f_a}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_y}}+\frac{2v_z}{w_z}\left[-(r_x+t_x)+\frac{w_x}{w_z}(r_z-t_z)\right]\label{eqn:df_aq_yq_y}\\
\frac{\partial^2{f_b}}{\partial{q_y^2}}&=\frac{-2}{w_z}\frac{\partial{f_b}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_y}}+\frac{2v_z}{w_z}\left[(r_y+t_y)+\frac{w_y}{w_z}(r_z-t_z)\right]\label{eqn:df_bq_yq_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_z^2}}&=\frac{-2}{w_z}\frac{\partial{f_a}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{2v_z}{w_z}\left[-(r_x-t_x)-\frac{w_x}{w_z}(r_z-t_z)\right]\label{eqn:df_aq_zq_z}\\
\frac{\partial^2{f_b}}{\partial{q_z^2}}&=\frac{-2}{w_z}\frac{\partial{f_b}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{2v_z}{w_z}\left[-(r_y-t_y)-\frac{w_y}{w_z}(r_z-t_z)\right]\label{eqn:df_bq_zq_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_w}\partial{q_x}}&=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_x}}+\frac{\partial{f_a}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_w}}\right]-\frac{2v_z}{w_z}\frac{w_x}{w_z}(r_y-t_y)\label{eqn:df_aq_wq_x}\\
\frac{\partial^2{f_b}}{\partial{q_w}\partial{q_x}}&=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_x}}+\frac{\partial{f_b}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left[-(r_z-t_z)-\frac{w_y}{w_z}(r_y-t_y)\right]\label{eqn:df_bq_wq_x}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_w}\partial{q_y}}&=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_y}}+\frac{\partial{f_a}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left[(r_z-t_z)+\frac{w_x}{w_z}(r_x-t_x)\right]\label{eqn:df_aq_wq_y}\\
\frac{\partial^2{f_b}}{\partial{q_w}\partial{q_y}}&=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_y}}+\frac{\partial{f_b}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\frac{w_y}{w_z}(r_x-t_x)\label{eqn:df_bq_wq_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_w}\partial{q_z}}&=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{\partial{f_a}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_w}}\right]-\frac{2v_z}{w_z}(r_y-t_y)\label{eqn:df_aq_wq_z}\\
\frac{\partial^2{f_b}}{\partial{q_w}\partial{q_z}}&=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{\partial{f_b}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}(r_x-t_x)\label{eqn:df_bq_wq_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_x}\partial{q_y}}&=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_y}}+\frac{\partial{f_a}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}(r_y-t_y)\label{eqn:df_aq_xq_y}\\
\frac{\partial^2{f_b}}{\partial{q_x}\partial{q_y}}&=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_y}}+\frac{\partial{f_b}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}(r_x-t_x)\label{eqn:df_bq_xq_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_x}\partial{q_z}}&=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{\partial{f_a}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}\left[(r_z-t_z)-\frac{w_x}{w_z}(r_x-t_x)\right]\label{eqn:df_aq_xq_z}\\
\frac{\partial^2{f_b}}{\partial{q_x}\partial{q_z}}&=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{\partial{f_b}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_x}}\right]-\frac{2v_z}{w_z}\frac{w_y}{w_z}(r_x-t_x)\label{eqn:df_bq_xq_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_y}\partial{q_z}}&=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{\partial{f_a}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_y}}\right]-\frac{2v_z}{w_z}\frac{w_x}{w_z}(r_y-t_y)\label{eqn:df_aq_yq_z}\\
\frac{\partial^2{f_b}}{\partial{q_y}\partial{q_z}}&=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{q_z}}+\frac{\partial{f_b}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{q_y}}\right]+\frac{2v_z}{w_z}\left[(r_z-t_z)-\frac{w_y}{w_z}(r_y-t_y)\right]\label{eqn:df_bq_yq_z}
\end{align}
When $q=r_x,r_y,r_z$ and the other variable is $p=q_w,q_x,q_y,q_z$ then,
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_w}\partial{r_x}}&=-\frac{\partial^2{f_a}}{\partial{q_w}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left(q_0+\frac{w_x}{w_z}q_2\right)\label{eqn:df_aq_wr_x}\\
\frac{\partial^2{f_b}}{\partial{q_w}\partial{r_x}}&=-\frac{\partial^2{f_b}}{\partial{q_w}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left(q_3+\frac{w_x}{w_z}q_2\right)\label{eqn:df_bq_wr_x}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_w}\partial{r_y}}&=-\frac{\partial^2{f_a}}{\partial{q_w}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left(-q_3-\frac{w_x}{w_z}q_1\right)\label{eqn:df_aq_wr_y}\\
\frac{\partial^2{f_b}}{\partial{q_w}\partial{r_y}}&=-\frac{\partial^2{f_b}}{\partial{q_w}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left(q_0-\frac{w_x}{w_z}q_1\right)\label{eqn:df_bq_wr_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_w}\partial{r_z}}&=-\frac{\partial^2{f_a}}{\partial{q_w}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left(q_2-\frac{w_x}{w_z}q_0\right)\label{eqn:df_aq_wr_z}\\
\frac{\partial^2{f_b}}{\partial{q_w}\partial{r_z}}&=-\frac{\partial^2{f_b}}{\partial{q_w}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_w}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_w}}\right]+\frac{2v_z}{w_z}\left(-q_1-\frac{w_x}{w_z}q_0\right)\label{eqn:df_bq_wr_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_x}\partial{r_x}}&=-\frac{\partial^2{f_a}}{\partial{q_x}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}\left(q_1-\frac{w_x}{w_z}q_3\right)\label{eqn:df_aq_xr_x}\\
\frac{\partial^2{f_b}}{\partial{q_x}\partial{r_x}}&=-\frac{\partial^2{f_b}}{\partial{q_x}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_b}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}\left(q_2-\frac{w_x}{w_z}q_3\right)\label{eqn:df_bq_xr_x}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_x}\partial{r_y}}&=-\frac{\partial^2{f_a}}{\partial{q_x}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}\left(q_2-\frac{w_x}{w_z}q_0\right)\label{eqn:df_aq_xr_y}\\
\frac{\partial^2{f_b}}{\partial{q_x}\partial{r_y}}&=-\frac{\partial^2{f_b}}{\partial{q_x}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_b}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}\left(-q_1-\frac{w_x}{w_z}q_0\right)\label{eqn:df_bq_xr_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_x}\partial{r_z}}&=-\frac{\partial^2{f_a}}{\partial{q_x}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}\left(q_3+\frac{w_x}{w_z}q_1\right)\label{eqn:df_aq_xr_z}\\
\frac{\partial^2{f_b}}{\partial{q_x}\partial{r_z}}&=-\frac{\partial^2{f_b}}{\partial{q_x}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_x}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_b}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_x}}\right]+\frac{2v_z}{w_z}\left(-q_0+\frac{w_x}{w_z}q_1\right)\label{eqn:df_bq_xr_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_y}\partial{r_x}}&=-\frac{\partial^2{f_a}}{\partial{q_y}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_y}}\right]+\frac{2v_z}{w_z}\left(-q_2+\frac{w_x}{w_z}q_0\right)\label{eqn:df_aq_yr_x}\\
\frac{\partial^2{f_b}}{\partial{q_y}\partial{r_x}}&=-\frac{\partial^2{f_b}}{\partial{q_y}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_b}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_y}}\right]+\frac{2v_z}{w_z}\left(q_1+\frac{w_x}{w_z}q_0\right)\label{eqn:df_bq_yr_x}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_y}\partial{r_y}}&=-\frac{\partial^2{f_a}}{\partial{q_y}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_y}}\right]+\frac{2v_z}{w_z}\left(q_1-\frac{w_x}{w_z}q_3\right)\label{eqn:df_aq_yr_y}\\
\frac{\partial^2{f_b}}{\partial{q_y}\partial{r_y}}&=-\frac{\partial^2{f_b}}{\partial{q_y}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_b}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_y}}\right]+\frac{2v_z}{w_z}\left(q_2-\frac{w_x}{w_z}q_3\right)\label{eqn:df_bq_yr_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_y}\partial{r_z}}&=-\frac{\partial^2{f_a}}{\partial{q_y}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_y}}\right]+\frac{2v_z}{w_z}\left(q_0+\frac{w_x}{w_z}q_2\right)\label{eqn:df_aq_yr_z}\\
\frac{\partial^2{f_b}}{\partial{q_y}\partial{r_z}}&=-\frac{\partial^2{f_b}}{\partial{q_y}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_y}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_b}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_y}}\right]+\frac{2v_z}{w_z}\left(q_3+\frac{w_x}{w_z}q_2\right)\label{eqn:df_bq_yr_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_z}\partial{r_x}}&=-\frac{\partial^2{f_a}}{\partial{q_z}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_z}}\right]+\frac{2v_z}{w_z}\left(-q_3-\frac{w_x}{w_z}q_1\right)\label{eqn:df_aq_zr_x}\\
\frac{\partial^2{f_b}}{\partial{q_z}\partial{r_x}}&=-\frac{\partial^2{f_b}}{\partial{q_z}\partial{t_x}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{r_x}}+\frac{\partial{f_b}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{q_z}}\right]+\frac{2v_z}{w_z}\left(q_0-\frac{w_x}{w_z}q_1\right)\label{eqn:df_bq_zr_x}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_z}\partial{r_y}}&=-\frac{\partial^2{f_a}}{\partial{q_z}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_z}}\right]+\frac{2v_z}{w_z}\left(-q_0-\frac{w_x}{w_z}q_2\right)\label{eqn:df_aq_zr_y}\\
\frac{\partial^2{f_b}}{\partial{q_z}\partial{r_y}}&=-\frac{\partial^2{f_b}}{\partial{q_z}\partial{t_y}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_b}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{q_z}}\right]+\frac{2v_z}{w_z}\left(-q_3-\frac{w_x}{w_z}q_2\right)\label{eqn:df_bq_zr_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{q_z}\partial{r_z}}&=-\frac{\partial^2{f_a}}{\partial{q_z}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_a}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_z}}\right]+\frac{2v_z}{w_z}\left(q_1-\frac{w_x}{w_z}q_3\right)\label{eqn:df_aq_zr_z}\\
\frac{\partial^2{f_b}}{\partial{q_z}\partial{r_z}}&=-\frac{\partial^2{f_b}}{\partial{q_z}\partial{t_z}}=\frac{-1}{w_z}\left[\frac{\partial{f_b}}{\partial{q_z}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_b}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{q_z}}\right]+\frac{2v_z}{w_z}\left(q_2-\frac{w_x}{w_z}q_3\right)\label{eqn:df_bq_zr_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{r_x}^2}&=\frac{\partial^2{f_a}}{\partial{t_x}^2}=-\frac{\partial^2{f_a}}{\partial{r_x}\partial{t_x}}=\frac{-2}{w_z}\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{r_x}}\label{eqn:df_ar_xr_x}\\
\frac{\partial^2{f_b}}{\partial{r_x}^2}&=\frac{\partial^2{f_b}}{\partial{t_x}^2}=-\frac{\partial^2{f_b}}{\partial{r_x}\partial{t_x}}=\frac{-2}{w_z}\frac{\partial{f_b}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{r_x}}\label{eqn:df_br_xr_x}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{r_y}^2}&=\frac{\partial^2{f_a}}{\partial{t_y}^2}=-\frac{\partial^2{f_a}}{\partial{r_y}\partial{t_y}}=\frac{-2}{w_z}\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{r_y}}\label{eqn:df_ar_yr_y}\\
\frac{\partial^2{f_b}}{\partial{r_y}^2}&=\frac{\partial^2{f_b}}{\partial{t_y}^2}=-\frac{\partial^2{f_b}}{\partial{r_y}\partial{t_y}}=\frac{-2}{w_z}\frac{\partial{f_b}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{r_y}}\label{eqn:df_br_yr_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{r_z}^2}&=\frac{\partial^2{f_a}}{\partial{t_z}^2}=-\frac{\partial^2{f_a}}{\partial{r_z}\partial{t_z}}=\frac{-2}{w_z}\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{r_z}}\label{eqn:df_ar_zr_z}\\
\frac{\partial^2{f_b}}{\partial{r_z}^2}&=\frac{\partial^2{f_b}}{\partial{t_z}^2}=-\frac{\partial^2{f_b}}{\partial{r_z}\partial{t_z}}=\frac{-2}{w_z}\frac{\partial{f_b}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{r_z}}\label{eqn:df_br_zr_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{r_x}\partial{r_y}}&=\frac{\partial^2{f_a}}{\partial{t_x}\partial{t_y}}=-\frac{\partial^2{f_a}}{\partial{r_x}\partial{t_y}}=-\frac{\partial^2{f_a}}{\partial{t_x}\partial{r_y}}=\frac{-1}{w_z}\left(\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{r_x}}\right)\label{eqn:df_ar_xr_y}\\
\frac{\partial^2{f_b}}{\partial{r_x}\partial{r_y}}&=\frac{\partial^2{f_b}}{\partial{t_x}\partial{t_y}}=-\frac{\partial^2{f_b}}{\partial{r_x}\partial{t_y}}=-\frac{\partial^2{f_b}}{\partial{t_x}\partial{r_y}}=\frac{-1}{w_z}\left(\frac{\partial{f_b}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{r_y}}+\frac{\partial{f_b}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{r_x}}\right)\label{eqn:df_br_xr_y}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{r_x}\partial{r_z}}&=\frac{\partial^2{f_a}}{\partial{t_x}\partial{t_z}}=-\frac{\partial^2{f_a}}{\partial{r_x}\partial{t_z}}=-\frac{\partial^2{f_a}}{\partial{t_x}\partial{r_z}}=\frac{-1}{w_z}\left(\frac{\partial{f_a}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{r_x}}\right)\label{eqn:df_ar_xr_z}\\
\frac{\partial^2{f_b}}{\partial{r_x}\partial{r_z}}&=\frac{\partial^2{f_b}}{\partial{t_x}\partial{t_z}}=-\frac{\partial^2{f_b}}{\partial{r_x}\partial{t_z}}=-\frac{\partial^2{f_b}}{\partial{t_x}\partial{r_z}}=\frac{-1}{w_z}\left(\frac{\partial{f_b}}{\partial{r_x}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_b}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{r_x}}\right)\label{eqn:df_br_xr_z}
\end{align}
\begin{align}
\frac{\partial^2{f_a}}{\partial{r_y}\partial{r_z}}&=\frac{\partial^2{f_a}}{\partial{t_y}\partial{t_z}}=-\frac{\partial^2{f_a}}{\partial{r_y}\partial{t_z}}=-\frac{\partial^2{f_a}}{\partial{t_y}\partial{r_z}}=\frac{-1}{w_z}\left(\frac{\partial{f_a}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_a}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{r_y}}\right)\label{eqn:df_ar_yr_z}\\
\frac{\partial^2{f_b}}{\partial{r_y}\partial{r_z}}&=\frac{\partial^2{f_b}}{\partial{t_y}\partial{t_z}}=-\frac{\partial^2{f_b}}{\partial{r_y}\partial{t_z}}=-\frac{\partial^2{f_b}}{\partial{t_y}\partial{r_z}}=\frac{-1}{w_z}\left(\frac{\partial{f_b}}{\partial{r_y}}\frac{\partial{w_z}}{\partial{r_z}}+\frac{\partial{f_b}}{\partial{r_z}}\frac{\partial{w_z}}{\partial{r_y}}\right)\label{eqn:df_br_yr_z}
\end{align}
\section{Formulation of Global Mathematical Model}
Up to this point all the equations derived have assumed only one object point, one image point with a single image and camera.  In this simple case there are (at most) fourteen unknown variables and two independent equations.  In order to continue there either needs to be assumptions made about some of the variables or more independent equations need to be found, or both.  In most cases, there are many image points with more than one image with correlating object points. Generally it is assumed that the image points in each image are known and it is desired to find the locations of the object points along with the camera and image parameters. Of course this process can be reversed and given the object points the image points can be found.  

Each object point with a correlating image point and camera generates two independent functons.  If an other image is taken of the same object point then an additonal set of independent equations are formed.  The point to be made is that there will be many independent functions, all of which need to be satisfied simultaneously.  Before they can be satisfied they must be combined and weighted relative to each other into one singe equation.  When this single global error function is minimized the best overall solution to all the indvidual functions will be found.

\subsection{Global Collinearity Equations}
Assume that every object point correlates to an image point for each image, or in other words every object point is visible by each image. Also, let there be $M$ images and $N$ object points with $K$ cameras. Each of the individual (local) collinearity equations with their respective indpendent variables are enumerated with the following indices,
\begin{itemize}
\item $\mathsf{m}$ = image index
\item $\mathsf{n}$ = object index
\item $\mathsf{k}$ = camera index
\end{itemize} 
How each variable is associated with the indices and the total number of variables is given in table \ref{tab:DescriptionOfIndependentVariableIndices}. 
\begin{table}
	\caption{Description of Independent Variable Indices}
	\centering
	\begin{tabular}{|l|l|l|l|}
	\hline
	Variables  & Index & No. Variables & Description\\
	\hline
	$t_{x\mathsf{m}}$, $t_{y\mathsf{m}}$, $t_{z\mathsf{m}}$,  & $\mathsf{m}$ & 6M & image index\\ 
	$\alpha_{x\mathsf{m}}$, $\alpha_{y\mathsf{m}}$, $\alpha_{z\mathsf{m}}$,& & & \\ 
	\hline
	$r_{x\mathsf{n}}$, $r_{y\mathsf{n}}$, $r_{z\mathsf{n}}$ & $\mathsf{n}$ & 3N & object index \\
	\hline
	$v_{x\mathsf{k_m}}$, $v_{y\mathsf{k_m}}$, $v_{z\mathsf{k_m}}$ & $\mathsf{k_m}$ & 3K & camera index \\
	\hline
	$u_{x\mathsf{m,n}}$, $u_{y\mathsf{m,n}}$ & $\mathsf{m,n}$ & 2MN & object/image index \\
	\hline
	\end{tabular}
	\label{tab:DescriptionOfIndependentVariableIndices}
\end{table}
The camera index $\mathsf{k_m}$ is dependent on the image index $\mathsf{m}$ since each image is taken with a single camera. The variable $\mathsf{k_m}$ can be implemented as a lookup table that tells which image was taken with what camera.  Note that the inverse of this question is not unique, since there can be many images taken with the same camera.

All of the individual functions are combined with an RMS type combining to get the following global error function,
\begin{align}
\varepsilon(q) = \sqrt{\frac{1}{MN}\sum_{\mathsf{m,n}}{f(p_{\mathsf{m,n}})^2}}\label{eqn:errorfcn}
\end{align}
where $q$ is a global set of independent variables and $p_\mathsf{m,n}$ is a local set of independent variables of each of the individual collinearity equations. The global set $q$ is the \emph{union} of all the local variables $p_\mathsf{m,n}$ or $q=\bigcup{p_\mathsf{m,n}}$ where,
\begin{align}
p_\mathsf{m,n}=\{r_{x\mathsf{n}},r_{y\mathsf{n}},r_{z\mathsf{n}},t_{x\mathsf{m}},t_{y\mathsf{m}},t_{z\mathsf{m}},v_{x\mathsf{k_m}},v_{y\mathsf{k_m}},v_{z\mathsf{k_m}},\alpha_{x\mathsf{m}},\alpha_{y\mathsf{m}},\alpha_{z\mathsf{m}},u_{x\mathsf{m,n}},u_{y\mathsf{m,n}}\}\label{eqn:photogramset_sub}
\end{align}
The reason that it is a \emph{union} of the local variables is that some variables in one local set will be the same variables as in another set. These variables need to be counted as a single variable in the global problem.  It is this fact that makes it possible to solve the photogrammetry problem, if there is no reduction in the number of variables in this manner there would always be more unknowns than equations.  The way this is actually accomplished or implemented is by the use of the three separate indices $\mathsf{m,n,k}$ to enumerate the variables, which is dependent on the assumption that every object point is visible by each image.  It should be emphasized that how the local variables correlated to the global variables is assumed known.  This is not an easy problem to solve, it can be solved manually by inspection or with other various algorithms such as stereo matchers and that use statistics and probability.

Once all of the local variables have been unified into a global set, they need to be separated into two subsets of known $q_k$ and unknown $q_u$ variables.  Only the \emph{subset} of unknown variables will be used in the optimization routines since obviously these are the only variables that need to be solved for. For example, if the image parameters $\vec{u}$ are assumed known then the subset of unknown $q_u$ is found to be,
\begin{align}
q_u&=\bigcup_{\mathsf{m,n,k_m}}\{r_{x\mathsf{n}},r_{y\mathsf{n}},r_{z\mathsf{n}},t_{x\mathsf{m}},t_{y\mathsf{m}},t_{z\mathsf{m}},v_{x\mathsf{k_m}},v_{y\mathsf{k_m}},v_{z\mathsf{k_m}},\alpha_{x\mathsf{m}},\alpha_{y\mathsf{m}},\alpha_{z\mathsf{m}}\}\\
\begin{split}
&=\{r_{x\mathsf{1}}\ldots{}r_{x\mathsf{N}},r_{y\mathsf{1}}\ldots{}r_{y\mathsf{N}},r_{z\mathsf{1}}\ldots{}r_{y\mathsf{N}},t_{x\mathsf{1}}\ldots{}t_{x\mathsf{M}},t_{y\mathsf{1}}\ldots{}t_{y\mathsf{M}},t_{z\mathsf{1}}\ldots{}t_{y\mathsf{M}},\\ &\qquad\alpha_{x\mathsf{1}}\ldots{}\alpha_{x\mathsf{M}},\alpha_{y\mathsf{1}}\ldots{}\alpha_{y\mathsf{M}},\alpha_{z\mathsf{1}}\ldots{}\alpha_{z\mathsf{M}},v_{x\mathsf{1}}\ldots{}v_{x\mathsf{K}},v_{y\mathsf{1}}\ldots{}v_{y\mathsf{K}},v_{z\mathsf{1}}\ldots{}v_{z\mathsf{K}}
\}
\end{split}
\end{align}
From the above discussion, it can be concluded that there are two basic ways of reducing the number of variables compared to the number of equations.  One is by combining the local variables into a global subset where variables in one local set will be the same as in another set.  The other is by assuming the variables are known, which will have to be obtained by some outside means or some other method. 

The local collinearity equation of (\ref{eqn:collinearcombine}) is enumerated with the subscripts as,
\begin{align}
f(p_{\mathsf{m},\mathsf{n}})=\sqrt{\frac{f_a(p_{\mathsf{m},\mathsf{n}})^2+f_b(p_{\mathsf{m},\mathsf{n}})^2}{2}}\label{eqn:functionmn}
\end{align}
where,
\begin{align}
f_a(p_{\mathsf{m},\mathsf{n}}) &= u_{x\mathsf{m},\mathsf{n}}-v_{x\mathsf{k_m}}+v_{z\mathsf{k_m}}\frac{w_{x}(p_{\mathsf{m},\mathsf{n}})}{w_{z}(p_{\mathsf{m},\mathsf{n}})}\label{eqn:collineara_sub}\\
f_b(p_{\mathsf{m},\mathsf{n}}) &= u_{y\mathsf{m},\mathsf{n}}-v_{y\mathsf{k_m}}+v_{z\mathsf{k_m}}\frac{w_{y}(p_{\mathsf{m},\mathsf{n}})}{w_{z}(p_{\mathsf{m},\mathsf{n}})}\label{eqn:collinearb_sub}
\end{align}
and where,
\begin{align}
\left[ \begin{array}{c}
w_{x}(p_{\mathsf{m},\mathsf{n}})\\w_{y}(p_{\mathsf{m},\mathsf{n}})\\w_{z}(p_{\mathsf{m},\mathsf{n}})
\end{array} \right]
=\left[ \begin{array}{ccc}
m_{11\mathsf{m}}(r_{x\mathsf{n}}-t_{x\mathsf{m}})+m_{12\mathsf{m}}(r_{y\mathsf{n}}-t_{y\mathsf{m}})+m_{13\mathsf{m}}(r_{z\mathsf{n}}-t_{z\mathsf{m}})\\
m_{21\mathsf{m}}(r_{x\mathsf{n}}-t_{x\mathsf{m}})+m_{22\mathsf{m}}(r_{y\mathsf{n}}-t_{y\mathsf{m}})+m_{23\mathsf{m}}(r_{z\mathsf{n}}-t_{z\mathsf{m}})\\
m_{31\mathsf{m}}(r_{x\mathsf{n}}-t_{x\mathsf{m}})+m_{32\mathsf{m}}(r_{y\mathsf{n}}-t_{y\mathsf{m}})+m_{33\mathsf{m}}(r_{z\mathsf{n}}-t_{z\mathsf{m}})
\end{array} \right]
\end{align}
Substituting (\ref{eqn:functionmn}) into (\ref{eqn:errorfcn}) yields,
\begin{align}
\varepsilon(q) = \sqrt{\frac{1}{2MN}\sum_{\mathsf{m,n}}[f_a(p_{\mathsf{m},\mathsf{n}})^2+f_b(p_{\mathsf{m},\mathsf{n}})^2}]\label{eqn:errorfcn2}
\end{align}
The gradient is calculated to be,
\begin{align}
\begin{split}	\nabla{\varepsilon(q)}&=
\sum_{\mathsf{n}}\left[\frac{\partial{\varepsilon(q)}}{\partial{r_{x\mathsf{n}}}}\hat{r}_{x\mathsf{n}}+\frac{\partial{\varepsilon(q)}}{\partial{r_{y\mathsf{n}}}}\hat{r}_{y\mathsf{n}}+\frac{\partial{\varepsilon(q)}}{\partial{r_{z\mathsf{n}}}}\hat{r}_{z\mathsf{n}}\right]\\
&+\sum_{\mathsf{m}}\left[\frac{\partial{\varepsilon(q)}}{\partial{t_{x\mathsf{m}}}}\hat{t}_{x\mathsf{m}}+\frac{\partial{\varepsilon(q)}}{\partial{t_{y\mathsf{m}}}}\hat{t}_{y\mathsf{m}}+\frac{\partial{\varepsilon(q)}}{\partial{t_{z\mathsf{m}}}}\hat{t}_{z\mathsf{m}}\right]\\
&+\sum_{\mathsf{m}}\left[\frac{\partial{\varepsilon(q)}}{\partial{\alpha_{x\mathsf{m}}}}\hat{\alpha}_{x\mathsf{m}}+\frac{\partial{\varepsilon(q)}}{\partial{\alpha_{y\mathsf{m}}}}\hat{\alpha}_{y\mathsf{m}}+\frac{\partial{\varepsilon(q)}}{\partial{\alpha_{z\mathsf{m}}}}\hat{\alpha}_{z\mathsf{m}}\right]\\
&+\sum_{\mathsf{k_m}}\left[\frac{\partial{\varepsilon(q)}}{\partial{v_{x\mathsf{k_m}}}}\hat{v}_{x\mathsf{k_m}}+\frac{\partial{\varepsilon(q)}}{\partial{v_{y\mathsf{k_m}}}}\hat{v}_{y\mathsf{k_m}}+\frac{\partial{\varepsilon(q)}}{\partial{v_{z\mathsf{k_m}}}}\hat{v}_{z\mathsf{k_m}}\right]\\
\end{split}
\end{align}
The derivative of (\ref{eqn:errorfcn}) is found to be,
\begin{align}
\varepsilon'(q)&= \frac{1}{MN\varepsilon(q)}\sum_{\mathsf{m},\mathsf{n}}f(p_\mathsf{m,n})f'(p_\mathsf{m,n})\label{eqn:partial_dg_p}
\end{align}
where,
\begin{align}
f'(p_\mathsf{m,n})&=\frac{f_a(p_\mathsf{m,n})f_a'(p_\mathsf{m,n})+f_b(p_\mathsf{m,n})f_b'(p_\mathsf{m,n})}{2f(p_\mathsf{m,n})}\label{eqn:partial_df_p}
\end{align}
Substituting, (\ref{eqn:partial_df_p}) into (\ref{eqn:partial_dg_p}) yields,
\begin{align}
\varepsilon'(q) = \frac{1}{2MN\varepsilon(q)}\sum_\mathsf{m,n}\left[f_a(p_\mathsf{m,n})f'_a(p_\mathsf{m,n})+f_b(p_\mathsf{m,n})f'_b(p_\mathsf{m,n})\right]\label{eqn:derrorfcn}
\end{align}

\subsection{Uniqueness of Solution}
Because of the assumption that every object point correlates to an image point for each image, the maximum number of variables $T$ in the problem is calculated to be,
\begin{align}
T = 3N+6M+3K+2MN
\end{align}
Under these assumptions, a simple counting exercise will show how many equations will be needed to solve for a given number of unknown variables. First, the number of unknown variables $U$ will be solved.  For every camera used, the unknown variables $v_{x\mathsf{k}}$,  $v_{y\mathsf{k}}$, and $v_{z\mathsf{k}}$ will be added, which will account for $3K$ unknowns.  Every image used will introduce the variables $t_{x\mathsf{m}}$,  $t_{y\mathsf{m}}$, and $t_{z\mathsf{m}}$, $\alpha_{x\mathsf{m}}$,  $\alpha_{y\mathsf{m}}$, and $\alpha_{z\mathsf{m}}$ which will add an additional $6M$ unknowns. Finally, for every object point the variables  $r_{x\mathsf{n}}$,  $r_{y\mathsf{n}}$, and $r_{z\mathsf{n}}$ is included and will add an additional $3N$ unknowns.  Of course, if any one of these variables are assumed known, they can be simply subtracted from the total number of unknowns.  This will be accounted for with the variable $U0$.  For example, if one image was located at the origin of the global coordinate system, then this would reduce the number of knowns by $3$. Also if the same image was oriented along one of the global coordinate system axes, then an additional $3$ unknowns could be taken off.  A number of additional assumptions can be made. The total number of unknowns $U$ is then given by,
\begin{align}
U=3K+6M+3N-U0\label{eqn:totalnumberofunknowns}
\end{align}
Every time an object point is added that correlates to an image point and an image, two independent equations are generated.  It is important to note the assumption that every object point is visible by every image is assumed in this case.  Therefore, the total number of independent equations $E$ is then given by,
\begin{align}
E = 2MN\label{eqn:totalnumberofeqns}
\end{align}
In order to obtain a unique solution to the problem, there needs to be more independent equations than unknowns, or
\begin{align}
E\ge{U}\label{eqn:eqnsgreaterthanunknowns}
\end{align}
Using (\ref{eqn:totalnumberofunknowns}) and (\ref{eqn:totalnumberofeqns}), the minimum number of object points that will guarantee the condition of (\ref{eqn:eqnsgreaterthanunknowns}) is given by,
\begin{align}
N_{min} = \mathrm{roundup}\left(\frac{3K+6M-U0}{2M-3}\right)
\end{align}

\subsection{Distortion Models}
There are four principal sources of departure from collinearity.  These are radial lens distortion, decentric lens distortion, image plane unflatness, and in plane distortion.
\subsubsection{Radial Lens Distortion}
Radial lens distortion can be represented by a polynomial with a given number of terms.  This model is shown in (\ref{eqn:faRLD}).
\begin{align}
	f_{aRLD} &= K_1(r^2-1)u_x+K_2(r^4)u_x+K_3(r^6-1)u_x\label{eqn:faRLD}\\
	f_{bRLD} &= K_1(r^2-1)u_y+K_2(r^4)u_y+K_3(r^6-1)u_y\label{eqn:fbRLD}
\end{align}
where,
\begin{align}
r = \sqrt{(u_x-v_x)^2+(u_y-v_y)^2}\label{eqn:r_distortion}
\end{align}
The more terms used will lead to more accurate results but will also increase the number of unknowns that need to be solved.  The lowest number of terms that will still provide the desired accuracy is desired to maximize computational efficiency.
It is also desired to know the partial derivatives with each of the variables.  In order to aid this effort, the following derivatives are calculated. First with respect to $u_x$ and $v_x$,
\begin{align}
	\frac{\partial r^2}{\partial{u_x}} &= -\frac{\partial r^2}{\partial{v_x}}=2(u_x-v_x)\\
	\frac{\partial r^4}{\partial{u_x}} &= -\frac{\partial r^4}{\partial{v_x}}=4r^2(u_x-v_x)\\
	\frac{\partial r^6}{\partial{u_x}} &= -\frac{\partial r^6}{\partial{v_x}}=6r^4(u_x-v_x)
\end{align}
and then with respect to $u_y$ and $v_y$,
\begin{align}
	\frac{\partial r^2}{\partial{u_y}} &= -\frac{\partial r^2}{\partial{v_y}}=2(u_y-v_y)\\
	\frac{\partial r^4}{\partial{u_y}} &= -\frac{\partial r^4}{\partial{v_y}}=4r^2(u_y-v_y)\\
	\frac{\partial r^6}{\partial{u_y}} &= -\frac{\partial r^6}{\partial{v_y}}=6r^4(u_y-v_y)
\end{align}

The partial derivatives are then found to be,
\begin{align}
	\frac{\partial f_{aRLD}}{\partial{K_1}} &=(r^2-1)u_x\\
	\frac{\partial f_{aRLD}}{\partial{K_2}} &=(r^4-1)u_x\\
	\frac{\partial f_{aRLD}}{\partial{K_3}} &=(r^6-1)u_x
\end{align}

\begin{align}
	\frac{\partial f_{bRLD}}{\partial{K_1}} &=(r^2-1)u_y\\
	\frac{\partial f_{bRLD}}{\partial{K_2}} &=(r^4-1)u_y\\
	\frac{\partial f_{bRLD}}{\partial{K_3}} &=(r^6-1)u_y
\end{align}

Let,
\begin{align}
P_k = k_1+2k_2r^2+3k_3r^4
\end{align}

\begin{align}
	\frac{\partial f_{aRLD}}{\partial{u_x}} &=u_x\frac{\partial r^2}{\partial{u_x}}P_k+\frac{f_{aRLD}}{u_x}\\
	\frac{\partial f_{aRLD}}{\partial{u_y}} &=u_x\frac{\partial r^2}{\partial{u_y}}P_k\\
	\frac{\partial f_{aRLD}}{\partial{v_x}} &=u_x\frac{\partial{r^2}}{\partial{v_x}}P_k=-u_x\frac{\partial{r^2}}{\partial{u_x}}P_k=\frac{f_{aRLD}}{u_x}-\frac{\partial f_{aRLD}}{\partial{u_x}}\\
	\frac{\partial f_{aRLD}}{\partial{v_y}}&=u_x\frac{\partial{r^2}}{\partial{v_y}}P_k=-u_x\frac{\partial{r^2}}{\partial{u_y}}P_k=-\frac{\partial f_{aRLD}}{\partial{u_y}}
\end{align}

\begin{align}
	\frac{\partial f_{bRLD}}{\partial{u_x}} &=u_y\frac{\partial r^2}{\partial{u_x}}P_k\\
	\frac{\partial f_{bRLD}}{\partial{u_y}} &=u_y\frac{\partial r^2}{\partial{u_y}}P_k+\frac{f_{bRLD}}{u_y}\\
	\frac{\partial f_{bRLD}}{\partial{v_x}}&=u_y\frac{\partial{r^2}}{\partial{v_x}}P_k=-u_y\frac{\partial{r^2}}{\partial{u_x}}P_k=-\frac{\partial f_{bRLD}}{\partial{u_x}}\\
	\frac{\partial f_{bRLD}}{\partial{v_y}} &=u_y\frac{\partial{r^2}}{\partial{v_y}}P_k=-u_y\frac{\partial{r^2}}{\partial{u_y}}P_k=\frac{f_{bRLD}}{u_y}-\frac{\partial f_{aRLD}}{\partial{u_y}}
\end{align}

\subsubsection{Decentric Distortion}
Decentering line distortion is manifest when the lens is not aligned correctly to the optical axis.  The equations that approximate this distortion is given in (\ref{eqn:faDLD}) and (\ref{eqn:fbDLD}).
\begin{align}
f_{aDLD} = P_1(r^2+2u_x^2)+2P_2u_xu_y\label{eqn:faDLD}\\
f_{bDLD} = P_2(r^2+2u_y^2)+2P_1u_xu_y\label{eqn:fbDLD}
\end{align}
where $r$ is defined in (\ref{eqn:r_distortion}).  The partial derivatives are calculated to be,
\begin{align}
\frac{\partial f_{aDLD}}{\partial{P_1}} &=r^2+2u_x^2\\
\frac{\partial f_{aDLD}}{\partial{P_2}} &=2u_xu_y\\
\frac{\partial f_{aDLD}}{\partial{u_x}} &=P_1\left(\frac{\partial{r^2}}{\partial{u_x}}+4u_x\right)+2P_2u_y\\
\frac{\partial f_{aDLD}}{\partial{u_y}} &=P_1\frac{\partial{r^2}}{\partial{u_y}}+2P_2u_x\\
\frac{\partial f_{aDLD}}{\partial{v_x}} &=P_1\frac{\partial{r^2}}{\partial{v_x}}=-P_1\frac{\partial{r^2}}{\partial{u_x}}\\
\frac{\partial f_{aDLD}}{\partial{v_y}} &=P_1\frac{\partial{r^2}}{\partial{v_y}}=-P_1\frac{\partial{r^2}}{\partial{u_y}}
\end{align}

\begin{align}
\frac{\partial f_{bDLD}}{\partial{P_1}} &=2u_xu_y = \frac{\partial f_{aDLD}}{\partial{P_2}}\\
\frac{\partial f_{bDLD}}{\partial{P_2}} &=r^2+2u_y^2\\
\frac{\partial f_{bDLD}}{\partial{u_x}} &=P_2\frac{\partial{r^2}}{\partial{u_x}}+2P_1u_y\\
\frac{\partial f_{bDLD}}{\partial{u_y}} &=P_2\left(\frac{\partial{r^2}}{\partial{u_y}}+4u_y\right)+2P_1u_x\\
\frac{\partial f_{bDLD}}{\partial{v_x}} &=P_2\frac{\partial{r^2}}{\partial{v_x}}=-P_2\frac{\partial{r^2}}{\partial{u_x}}\\
\frac{\partial f_{bDLD}}{\partial{v_y}} &=P_2\frac{\partial{r^2}}{\partial{v_y}}=-P_2\frac{\partial{r^2}}{\partial{u_y}}
\end{align}

\subsubsection{Image Plane Unflatness}
The way to compensate for image plane unflatness is to measure the surface of image plane and then fit a polynomial to it.

\subsubsection{In-Plane Distortion}
In-plane distortions are usually manifest by differential scaling between x and image coordinates.  It can also cause the x and y axis to to be non-orthogonal.  These distortions are denoted as affine deformations and can be mathematically modeled according to (\ref{eqn:faAD}) and (\ref{eqn:fbAD}).  

\begin{align}
	f_{aAD} &= -A_1u_x+A_2u_y\label{eqn:faAD}\\
	f_{bAD} &= -A_1u_y\label{eqn:fbAD}
\end{align}

Where the partial derivatives are calculated to be,
\begin{align}
\frac{\partial f_{aAD}}{\partial{A_1}} &= -u_x\\
\frac{\partial f_{aAD}}{\partial{A_2}} &= u_y\\
\frac{\partial f_{aAD}}{\partial{u_x}} &= -A_1\\
\frac{\partial f_{aAD}}{\partial{u_y}} &= A_2
\end{align}

\begin{align}
\frac{\partial f_{bAD}}{\partial{A_1}} &= u_y\\
\frac{\partial f_{bAD}}{\partial{A_2}} &= 0\\
\frac{\partial f_{bAD}}{\partial{u_x}} &= 0\\
\frac{\partial f_{bAD}}{\partial{u_y}} &= A_1
\end{align}

\subsubsection{Extended Collinearity Equations}
The total distortions is found by superimposing each of the contributions of the distortions to get,
\begin{align}
f_{aD} &= f_{aRLD} + f_{aDLD} + f_{aAD} + f_{etc.}\label{eqn:faD}\\
f_{bD} &= f_{bRLD} + f_{bDLD} + f_{bAD} + f_{etc.}\label{eqn:fbD}\\
\end{align}

The distortions are included into the collinearity equations (\ref{eqn:collinear1}) and (\ref{eqn:collinear2}) by superposition to get,
\begin{align}
f_a(p) &= u_x-v_x+v_z\frac{m_{11}(r_x-t_x)+m_{12}(r_y-t_y)+m_{13}(r_z-t_z)}{m_{31}(r_x-t_x)+m_{32}(r_y-t_y)+m_{33}(r_z-t_z)}-f_{aD}\label{eqn:collinearD1}\\
f_b(p) &= u_y-v_y+v_z\frac{m_{21}(r_x-t_x)+m_{22}(r_y-t_y)+m_{23}(r_z-t_z)}{m_{31}(r_x-t_x)+m_{32}(r_y-t_y)+m_{33}(r_z-t_z)}-f_{bD}\label{eqn:collinearD2} 
\end{align}

\subsection{Solution of Global Collinearity Equations}
The photogrammetry problem can be solved by minimizing (\ref{eqn:errorfcn}}) and by making use of gradient information in (\ref{eqn:derrorfcn}).  There are different methods that can be used to solve this problem including,
\begin{itemize}
\item simplex
\item steepest descent
\item conjugate gradient
\item variable metric
\end{itemize} 

\section{Gradient in terms of non-linear residuals}
Suppose you have the following basis function or residual,
\begin{align}
	r_m(\vec{x})
\end{align}
where $m=1\ldots\infty$ and
\begin{align}
	\vec{x}\equiv\sum_{i}x_i\hat{x}_i
\end{align}
Now we form a general cost function $f$ by taking the p-norm,
\begin{align}
	f &= \left[\sum_m{r_m^p(\vec{x})}\right]^\frac{1}{p}\label{eqn:costfcn}
\end{align}
The gradient of a sclar funciton $f(\vec{x})$ is,
\begin{align}
	\nabla f = \sum_i{\frac{\partial f}{\partial x_i}}\hat{x_i}\label{eqn:gradientf}
\end{align}
Substituting (\ref{eqn:costfcn}) into (\ref{eqn:gradientf}) yields,
\begin{align}
	\nabla f &= \sum_i{\frac{\partial}{\partial x_i}}\left[\sum_m{r_m^p}\right]^\frac{1}{p}\hat{x_i}\\
	&= \sum_i\frac{1}{p\left[\sum_m{r_m^p}\right]^{\frac{-1}{p}}}\left[\sum_m{p\;r_m\frac{\partial r_m}{\partial x_i}}\right]\hat{x_i}\\
	&= \frac{1}{f}\sum_i\sum_m\left[{r_m\frac{\partial r_m}{\partial x_i}}\right]\hat{x_i}\\
	&= \frac{1}{f}\sum_m\sum_i\left[{r_m\frac{\partial r_m}{\partial x_i}}\right]\hat{x_i}
\end{align}



\end{document}